%dXt = sigma(t)dW 
%sigma(t)=OU
clear all
% index J for X
alphaK=20; betaK=15; 
% time T
T = 1/2;
% time mesh 
num=50; tspan = linspace(0,T,num+1);
%sigma
sig=sigmahandle; a=1;b=1;m=1;sig0=1;
s  =mihandle;
%options = odeset('AbsTol',1e-12,'RelTol',1e-6);
options = [];
row=(betaK+1)*(alphaK+1); X0=zeros(row,1); 
% solver of X
disp('solving X')
tic
[time,X] = ode45(@randucoudw,tspan,X0,options,T,a,b,m,sig0,alphaK,betaK,sig,s);

%varitest(time,X.*X,T,'uncorrelated','oudt',num)

hold on
%%%the WCE price
price=zeros(51,1);
Npath=10000;
eta=randn(Npath,betaK);
xi=randn(Npath, alphaK);

for iter=1:51
asset=X(iter,1)*ones(Npath,1);
for jj=1:betaK
    asset=asset+X(iter,jj+1)*eta(:,jj);
end
for ii=1:alphaK
    asset=asset+X(iter,ii*(betaK+1)+1)*xi(:,ii);
    for jj=1:betaK
        asset=asset+X(iter,ii*(betaK+1)+1+jj)*xi(:,ii).*eta(:,jj);
    end
end
price(iter)=sum(abs(asset))/Npath;
end
toc
plot(time, price);
%%% MC price
tic
mcpart=zeros(51,1);
for pn=1:Npath
    samplepath=getapath(5001,'sigmadW');
    sigmadW=abs(samplepath); 
    for jj=1:51
        mcpart(jj)=mcpart(jj)+sigmadW((jj-1)*100+1);
    end
end
toc
mcpart=mcpart/Npath;
plot(time,mcpart,'r')